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ABSTRACT 

The observation of blank fields, regions of the sky devoid of stars down to a given thresh- 
old magnitude, constitutes one of the typical important calibration procedures required for 
the proper reduction of astronomical data obtained in imaging mode. This work describes a 
method, based on the use of the Delaunay triangulation on the surface of a sphere, that allows 
the easy generation of blank fields catalogues. In addition to that, a new tool named TESELA, 
accessible through the WEB, has been created to facilitate the user to retrieve, and visualise 
using the VO-tool Aladin, the blank fields available near a given position in the sky. 

Key words: methods: data analysis - methods: numerical - Virtual observatory tools 
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1 INTRODUCTION 

Scientific research progress is critically based on the correct ex- 
ploitation of data obtained at the limit of the technological capa- 
bilities. Thus prior to data analysis and interpretation, the quality 
of the data treatment must be assured in order to guarantee that 
the information content in that data is preserved. In observational 
Astrophysics, and in particular imaging with ground-based tele- 
scopes, data treatment is performed following data reduction pro- 
cedures. The main goal of such reduction processes is to minimize 
the influence of data acquisition imperfections on the estimation of 
the desired astronomical measurements (see e.g.[Gilliland 1992]for 
a short review on noise sources and reduction processes of CCD 
data). For this purpose, one must perform appropriate manipula- 
tions of data and calibration frames, the latter properly planned and 
obtained to facilitate the reduction procedure. In this sense, image 
flatfielding and sky subtraction constitute two of the most common 
and important reduction steps. Their relevance should not be un- 
derestimated, since inadequate flatfielding or sky subtraction eas- 
ily leads to the introduction of systematic uncertainties in the data. 
Contrary to random errors, which can be controlled by applying 
normal statistical methods, the systematic uncertainties can propa- 
gate hidden within the arithmetically manipulated data to the final 
measurements, where their impact may considerably bias their in- 
terpretation. 

Concerning image flatfielding, the corresponding calibration 
frames needed to proceed through this reduction step are typi- 
cally divided into two categories: images intended to provide high- 
frequency scale sensitivity variations (the pixel-to-pixel response), 
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and images whose goal is to facilitate the removal of the two- 
dimensional low-frequency scale sensitivity variations of the detec- 
tors. The former are usually obtained from continuum lamps (either 
within the instrument itself or by illuminating a flat screen or the 
inner dome of the telescope), although the colour of such lamps 
are not expected to match the spectral energy distributions of the 
science targets, and ideally dark night sky flats, obtained by point- 
ing the telescope to a blank field should be observed instead. Since 
this is quite expensive in terms of observing time, this approach is 
very rarely used. The low-frequency flatfields are usually twilight 
exposures. The problem here is that the time during which the sky 
surface brightness is high enough to clearly overpass the signal of 
any star in the field of view, but not too high to saturate the frames, 
are placed on two time windows every day, one at the beginning of 
the evening twilight and the other at the end of the morning twi- 
light. The observational strategy typically consist on obtaining a 
series of images, shifting the telescope position a few arcseconds 
between the different exposures in order to avoid the bright stars in 
the field of view to appear in the same detector pixels. Those pixels 
are afterwards masked when combining the individual exposures 
during the reduction procedure. In any case very bright stars need 
to be avoided since the point spread function of these objects can 
introduce spikes and artifacts that are not so easily removed dur- 
ing the reduction of the flatfield images. For this purpose the use of 
blank fields is the best option. 

Sky subtraction cannot always be performed by measuring the 
night sky level in the same frame where the scientific targets are 
present. This is especially true when the target dimensions are com- 
parable to the detector field of view. In that circumstance separate 
sky frames are observed, by pointing the telescope to a position 
devoid of stars as much as possible. For the sky subtraction proce- 
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dure to be successful, it is essential that these separate sky images 
are obtained under identical observing conditions, i.e. very close in 
time and sky position, since this is the only way to prevent varying 
observing conditions that modify the sky surface brightness. 

From the above discussion it is clear that the observation of 
sky regions devoid of stars down to a given magnitude is a very 
important aspect to be taken into account when carrying out astro- 
nomical observations. To date, no systematic catalogue is available 
providing the location in the sky of blank fields for varying limit- 
ing magnitude^ With modern instrumentation intended to be used 
in large imaging surveys, the problem is expected to become more 
demanding. The situation is also starting to be important in spec- 
troscopic observations, since the advent of integral field units with 
high multiplexing capabilities (and thus, increasingly larger field of 
view), are expected to be one of the most common instrumentation 
in ground-based observatories. 

The work presented in this paper describes a method that helps 
to determine the availability of blank fields in any region of the 
celestial sphere, based on the use of the Delaunay triangulation. 
The method has been employed to generate catalogues of blank 
fields to varying limiting magnitudes. In order to facilitate the use 
of these catalogues, we have created TESELA, a new tool accessi- 
ble through the WEB which provides a simple interface that allows 
the user to retrieve the list of blank fields available near a given 
position in the sky. 

Section [2] of this paper describes the method and presents 
some statistical analysis of its application to the Tycho-2 (H0g et al. 
2000) catalogue. Section [3] describes the new tool TESELA. The 
detailed description of the tessellation of sky subregions is pre- 
sented in the Appendix. 



2 TESSELLATING THE SKY 



2.1 The Delaunay Triangulation 



The Delaunay triangulation (Delaunay 1934) is a subdivision of a 
geometric object (e.g. a surface or a volume) into a set of simplices. 
A simplex, or n- simplex to be specific, is the n-dimensional ana- 
logue of a triangle. More precisely, a simplex is the convex hull 
(convex envelope) of a set of (n + 1) points. 

In particular, for the Euclidean planar (2-dimensional) case, 
given a set of points, also called nodes, the Delaunay triangulation 
becomes a subdivision of the plane into triangles, whose vertices 
are nodes. For each of these triangles, it is possible to determine its 
associated circumcircle, the circle passing exactly through the three 
vertices of the triangle, and whose center, the circumcentre, can 
easily be computed as the intersection of the three perpendicular 
bisectors. Interestingly, in a Delaunay triangulation all the triangles 
satisfy the empty circumcircle interior property, which states that 
all the circumcircles are empty, i.e., there are no nodes inside any 
of the computed circumcircles. 



2.2 Applying the Delaunay triangulation to the celestial 
sphere 

Fortunately the Delaunay triangulation is not restricted to the Eu- 
clidean 2-dimensional case. It can be applied to other surfaces, in 
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particular to the 2-dimensional surface of a 3 -dimensional sphere. 
In this situation, the tessellation of the spherical surface is built 
with spherical triangles, i.e., triangles on the sphere whose sides 
are great circles. 

The empty circumcircle interior property of the Delaunay 
triangulation provides a straightforward method for a systematic 
search of regions in the celestial sphere free from stars. If one com- 
putes the Delaunay triangulation in the 2-dimensional surface of 
a sphere, using as nodes the location of the stars down to a given 
threshold visual magnitude (m), the above property guarantees that 
all the circumcircles are void of stars brighter than that magnitude 
(see Fig. [TJ. Thus, the circumdiameter of every circumcircle de- 
termines the maximum field of view that can be employed in that 
region of the sky as blank field. 

In this work we have applied the Delaunay triangulation to the 
Tycho-2 stellar catalogue. Tycho-2 contains astrometric and pho- 
tometric information for the 2.5 million brightest stars in the sky, 
and it is complete up to magnitude V = 11.5 mag. Photometric data 
consists of two pass-bands (Bt and Vr, close to Johnson B and V; 
Perryman et al. 1997). Typical uncertainties are 60 mas in posi- 
tion and 0. 1 mag in photometry. In addition to the main catalogue, 
we have used the first Tycho-2 supplement, which lists another 
17,588 bright stars from the Hipparcos and Tycho-1 Catalogues 
which are not in Tycho-2. 

In order to proceed with the triangulation, we have made use 
of STRIPACK ( |Renka|1997] >, a Fortran 77 software package that 
employs an incremental algorithm to build a Delaunay triangula- 
tion of a set of points on the surface of the unit sphere. For N 
nodes, the storage requirement for the triangulation is 13N integer 
storage locations in addition to 3N nodal coordinates. The compu- 
tation scales as O(NlogN). It is important to highlight that the 
original software was written using single-precision floating arith- 
metic, which turns out to be insufficient when dealing with astro- 
nomical coordinates with accuracy ~ 1 arcsec. For the work pre- 
sented here, we have modified the software code in order to use 
double -precision floating arithmetic, which guarantees the proper 
computation of the triangulation when working with star separa- 
tions approaching a few arc seconds. 

2.2.7 Tessellating a collection of smaller subregions 

In principle it is straightforward to obtain a list of blank fields for 
the whole celestial sphere by using as input a stellar catalogue in- 
cluding stars down to a given magnitude. However, considering that 
the number of stars grows very rapidly with increasing limiting 
magnitude, this approach may turn out to become nonviable, ei- 
ther in terms of computer memory or considering computing time. 
In Fig.[2]we show the time taken by our computer in function of the 
number of stars (nodes) to carry out the tessellation of the Tycho-2 
catalogue. Each dot corresponds to an increase of 0.5 mag in m y T , 
in the range from 6.0 to 10.0 mag. Although the CPU time depends 
on the characteristic of the machine used for the tessellation, the 
most important aspect is the exponential increase of the required 
time with mv T , mainly due to the exponential growth of the num- 
ber of stars (see Table [T}. In addition, since the tessellation of the 
whole celestial sphere demands memory storage (for stellar coordi- 
nates plus auxiliary variables) that scales with the number of stars, 
this approach requires excessively large memory capacity for faint 
magnitudes. 

In these circumstances our approach is to subdivide the whole 
celestial sphere in many smaller subregions (with some overlap 
among them), in order to individually compute the Delaunay tri- 
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Figure 1. Application of the Delaunay triangulation to the celestial sphere. Left: Sky region centered at the Orion constellation (a = 5 h 40 m and 5 = 0°, 
J2000.0), represented using a Lambert's equal area projection. All the stars brighter than my = 4.5 mag have been considered. Although each star is repre- 
sented with filled circles of size proportional to its brightness, the triangulation is carried out considering all the stars as nodes for the triangles, independently 
of their magnitude. A particular triangle is highlighted in red colour. Right:. Circumcircles corresponding to all the triangles derived in the previous triangula- 
tion. The red circle indicates the circumcircle associated to the previously highlighted triangle. Note that the empty circumcircle interior property guarantees 
that there are no stars inside any of the displayed circumcircles. 



25 



20 



15 



£ 10 



4.0 4.5 5.0 5.5 

log(N) 

Figure 2. CPU time consumed in function of the number of stars (nodes) 
used in the tessellation. Each point corresponds to an increasing my T from 
6 to 10 in steps of 0.5 mag. 



angulation in each of these smaller patches and, finally, merge the 
separate blank field lists. However the merging process is not im- 
mediate, since in this scenario the fact that the triangulation is com- 
puted in small subregions have an additional side effect that must be 
properly handled: when the subregion area is smaller than a single 
hemisphere, the convex envelope of the stars within that subregion 
has an interior zone (the subregion) and a larger exterior zone (the 
celestial sphere excluding the considered subregion). This means 
that the triangulation has boundary edges, i.e. triangles sides that 
are not shared with other triangles (see details in Renka 1997}. In 
practice, the problem arises when one computes the circumcircles 



of the boundary triangles, since in these cases those circumcircles 
may easily expand beyond the limits of the considered subregion, 
where no stars are present because they have not been included as 
input for the triangulation. A graphical illustration of this problem 
is shown in Fig. [3] 

Figure[3^ displays the initial triangulation obtained after com- 
puting the Delaunay triangulation of the subregion of the celes- 
tial sphere centered around the equatorial coordinates a — 6 h and 
S — 0° (J2000.0) within a radius of 30°, and using stars brighter 
than my = 4.5 mag (note that the Orion constellation is close to 
the center of the field, whereas the bright stars Sirius, Procyon and 
Aldebaran appear at the South-East, East and North- West, respec- 
tively). The sky region is represented using a Lambert's equal area 
projection ( Calabretta & Greisen 2002). The thick red line indicates 
the limit of the 30° radius. Once the triangulation has been com- 
puted, the next step is the calculation of the circumcircle associated 
to each triangle. Fig. [3]} shows all the circumcircles obtained for the 
triangulation derived in the previous step. The circumcircles that 
are fully circumscribed within the 30° radius are displayed with 
green colour, whereas the circumcircles that expand beyond that 
limit are represented in magenta. It is obvious that the latter will not 
be, in general, appropriate blank fields, since they have been com- 
puted assuming that no stars were present beyond the red thick line. 
However, even in these cases several possibilities can be envisaged 
in order to obtain valid blank fields starting from the information 
provided by the affected boundary triangles, while respecting the 
constraint imposed by the proximity of the boundary of the subre- 
gion. We discuss in Appendix |A| several of these possibilities. The 
result of following such approach is displayed in Fig.[5J;. 

Finally, the blank field lists obtained for the different sky sub- 
regions need to be merged. In our case, we have applied this pro- 
cedure when tessellating the sky for mv T = 10.5 and 11.0 mag. 
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Figure 3. Graphical illustration of the tessellation procedure when handling 
a subregion of the sky, in particular the portion of the celestial sphere (here 
represented using a Lambert's equal area projection) centered around the 
equatorial coordinates a = 6 h and 5 = 0° (J2000.0) within a radius of 30°, 
and using stars brighter than my = 4.5 mag. See explanation in the text 
(section |2TT) . 



In particular, we have divided the whole celestial sphere in 578 cir- 
cular subregions of 8° radius, where the centre of each subregion 
is separated by 6° in declination and by 6°/ cos 5 in right ascen- 
sion. Considering that the maximum blank field radius obtained for 
m v T — 10.0 mag is ~ 1° (see Table [TJ, this separation between 
subregions provides an excellent overlap among them. Within the 
overlapping areas there are duplicated blank fields and additional 
blank fields that were not permitted to expand beyond their cor- 
responding boundary limit. Consequently, to obtain the final cata- 
logue, we have removed both the repeated blank fields and those 
due to the border effect (typically inscribed within larger blank 
fields of neighbouring subregions). The latter were easy to identify 
because they circumscribe less than 3 stars. Thus, the final all-sky 
blank field catalogue is not affected by the number, size and overlap 
of the subregions that we have employed. 



2.2.2 Preparing the nodes for the triangulation 

As described above, Tycho-2 provides photometry data in two pass- 
bands (Bt and Vr). We have applied the Delaunay triangulation to 
each pass-band separately, and to the combination of both. To ob- 
tain the initial collections of stars, we have selected from both the 
main and the supplement Tycho-2 catalogues all stars with magni- 
tude lower than m in a given filter. In the case of using both filters, 
the star magnitude should be lower than m in either of the two fil- 
ters. 

These initial collections of stars are still not suitable to com- 
pute the Delaunay triangulation due to the presence of stars with 
insufficient separation. For that reason, we decided to "merge" into 
single objects all the stars closer than 1 arcsec, similar to a typ- 
ical seeing under good wheather conditions. The resulting visual 
magnitude for the combined objects was computed as the sum of 
the fluxes of the merged stars. The coordinates of the new objects 
were placed in the line connecting the merged stars, closer to the 
brightest star (using a weighting scheme dependent on the individ- 
ual brightness of the combined stars). This merging process does 
not reduce substantially the final number of stars, but it makes the 
triangulation process easier by removing unnecessary very small 
triangles for which the computations are prone to rounding errors. 



2.3 Results 

We have applied the triangulation method to the merged star lists as 
described above, with threshold magnitudes between 6 and 1 1 mag 
in steps of 0.5 mag. The resulting blank field catalogues are acces- 
sible through the TESELA tool (explained in section [3}. 

As an example, we present the results of the triangulation for 
Vt in Fig. [4] The different histograms correspond to the distribu- 
tions of the number of blank fields as a function of the blank field 
radius, for different threshold magnitudes mv T ■ The resulting his- 
tograms are positively skewed and although, as a first order ap- 
proach, they can approximately be fitted with a normal distribution 
(green curves), better fits are obtained using xl functions with v 
degrees of freedom (red curves). Note that the use of a chi-square 
law is just an empirical result, obtained after trying to fit different 
well-known skewed distributions, like lognormal and F U1)U2 distri- 
butions (the latter with v\ and v<i degrees of freedom), and finding 
that the best fits were obtained using F U1 ,^ 2 wh en v 2 — ^ oo, which 
is equivalent to use a xtt distribution (see e.g. Press et al.|2007 >. 
It is not the goal of this work to provide a physical justification for 
this behavior. 
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Figure 4. Histograms with the distributions of blank field radius, for different values of m y T , as indicated in the upper right inset of each diagram. Note that 
the x-axis range decreases as my T increases. The green curve is the fit to a normal distribution, whereas the red line indicates the corresponding fit to a xj 
distribution with v degrees of freedom, respectively. The coefficients of these parametric fits are given in Table [T] The best fit in all the diagrams is obtained 
for the xt distribution. 

Table 1. Results of the tessellation process using the Tycho-2 catalogue with the Vt filter. Down to my T = 10 mag, the Delaunay triangulation has been 
obtained using the resulting stellar catalogues corresponding to the whole celestial sphere. For my T = 10.5 and 11.0, the results correspond to combining 
the individual triangulation of smaller sky subregions, as explained in the text (see section [Z2.1) . The last five rows of this table indicate the coefficients 
corresponding to the linear regressions y = a + bx, where x = my T in all the cases, and y has been set to the logarithm (with base 10) of the parameters of 
the remaining columns where the linear fit makes sense. See description of the different columns in Sect.|2.3| 
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The quantitative description of the above results are presented 
in Table IT] The threshold magnitude my T is given in the first 
column. The number of stars iV* (i.e. nodes) and the number of 
blank field regions found Nbf are listed in the second and third 
columns, respectively. The fourth column indicates the median 
blank field radius pmedian and column (5) is the bin width Ah for 



the histogram distributions (some of them displayed in Fig.[4j em- 
ployed to derive the parametric fits given in the next columns. 
The amplitude A, mean radius p and standard deviation a p , ob- 
tained from the fit of the each histogram to a normal distribution 
of the form Aexp[— (p — p) 2 /(2af ) )] are listed in columns (6) 
to (8), respectively. The maximum blank field radius appears in 
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column (9). The coefficients of the fit to a chi- square distribution 
of the form B xt(p/f)> where B is the amplitude, / is a scaling 
factor for p, and v is the number of degrees of freedom, are given 
in columns (10) to (12). Finally, column (13) indicates the CPU 
time taken for the triangulation when tessellating the whole celes- 
tial sphere at once. 

As we have shown in Fig [2] the required CPU computing time 
increases exponentially with m. For this reason, we have derived 
the Delanuay triangulation for the whole celestial sphere at once 
just up to magnitude 10, and in overlaping smaller regions (as ex- 
plained in section |2!2T| ) for higher m. 



2.4 Analysis 

We have analysed the results of the triangulation for Vt- For Bt 
and the combination of both filters the results are very similar. 

Not surprisingly, there is an excellent correlation between 
most of the parameters listed in Table [T] and m y T . The correla- 
tions can be very well fitted using linear regressions of the form 
y = a-\-bx, being x = mv T and defining y as the base- 10 log- 
arithm of the considered parameter. The intercept a and slope b 
of the regressions (together with their associated uncertainties Aa 
and A b) are given at the bottom of each column. The final entry 
indicates the coefficient of determination r 2 , which in all the cases 
reveals the excellent correlation between the fitted data. 

It is also interesting to examine the characteristics of the de- 
rived blank fields at a given limiting magnitude. For example we 
have analyzed two illustrative diagrams corresponding to the case 
mv T = 11.0, which are represented in Fig. [5] Diagram, pja) dis- 
plays the variation in the number of regions as a function of the 
absolute value of galactic latitude /. Not surprisingly, the number 
of regions must decrease as the latitude increases, since the area 
of a spherical annulus (shown with the red symbols) is maximum 
at the galactic equator and tends to zero when \l\ approaches 90°. 
However, this variation is not enough to explain the difference be- 
tween the histogram and the red symbols, which clearly indicates 
that the number of regions is highly concentrated toward the galac- 
tic equator, which is simply the result of the higher stellar densitiy 
in that region. The diagram [5jb) shows that even though the num- 
ber of available blank fields unavoidably decreases with galactic 
latitude, their typical size, and the variation of sizes among them, 
increase with / approaching the galactic poles. 

The extrapolation to fainter magnitudes of the exponential 
variations of the parameters listed in Table [T] together with the un- 
avoidable increasing difficulty when approaching low galactic lat- 
itudes, strongly supports the need for an automatic tool that helps 
to identify suitable blank fields when observing with medium/large 
size telescopes. 



3 TESELA 

In order to provide an easy way to access the blank field catalogues, 
we have created TESELA, a WEB accesible tool, developed within 
the Spanish Virtual Observator)Q which can be publicly accessed 
through the following URL: 
|http ://sdc. cab .inta-csic.es/tesela| . 

The tool consists of a database containing the Tycho-2 stars 
and the already computed blank fields regions, and an user-friendly 
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5. Properties of the blank fields obtained 



Figure 5. Properties of the blank fields obtained in the case 
my T = 11.0 mag. Diagram (a): number of blank fields as a function of 
the absolute value of the galactic latitude. The red symbols show the rela- 
tive variation of the surface area on the celestial sphere with galactic lati- 
tude. Diagram (b): variation of the mean circumcircle radius as a function 
of the absolute value of the galactic latitude (the error bars indicate the root 
mean square exhibited by the radii at each galactic latitude bin). 



interface for accesing the data. TESELA allows the users to per- 
form a cone-search around a position in the sky, i.e. to obtain the list 
of blank fields available around a given sky position (right ascen- 
sion and declination) and within a fixed radius around that position. 
Through its search form (see Fig[6|, users can select the threshold 
magnitude m, the Tycho-2 filter to be used (Bt, Vt, or the com- 
bination of both), and define a minimum radius for the blank field 
regions. This last option may be especially important to fulfill ob- 
serving requirements, e.g., to ensure that the blank field is larger 
than the field of view of a given instrument. 

TESELA presents the result of the cone- search in two ta- 
bles, one for blank fields (RA, DEC, radius) and the other for the 
Tycho-2 stars (RA, DEC, jub t and/or mv T ) in the searching area. 
These tables can be downloaded in CSV format for further use. 
TESELA also provides users with the possibility of visu alizing the 
data. To do that, TESELA takes advantage of Aladirjj ( JBonnarel 
|et al.|2000| , a Virtual Observatory (VO) compliant software that al- 
lows users to visualise and analyse digitised astronomical images, 
and superimpose entries from astronomical catalogues or databases 
available from the VO services. Thanks to this connection with Al- 



|http: 



//svo . cab. inta-csic.es/ 



http: //aladin.u-strasbg. f r/ 



TESELA: blank fields for astronomical observations 7 




TESELA fif 



Sky Position (RA, DEC (J2000)): 
|12 19 3.0 +27 28 36.5 
Radius of field of view: 



TESELA: Search Form 



Examples: 

183.005 28.536 - Coordinates in decimal degrees 

12 19 3.0 +27 28 36.5 - Coordinates in RA(h,m,s) DEC(°,',") 



n radius of the blankfield: 



Threshold magnitude: 

Filter: 

| Tycho Bt&Vt 



(decimal degrees) 



(decimal degrees) 



Submit Query Reset Form 



Version 0.3 - Nov 2010 © CAB (INTA-CSIC) 



Home - Help Desk 



Figure 6. TESELA' s search form. Through it the user defines the value of 
several parameters used for the cone-search. This form is available at 

http : //sdc. cab. inta-csic .es/tesela 
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Figure 7. Example of visualization of the data with Aladin. The large red 
circle defines the searching area, Tycho-2 stars are shown by red diamonds, 
the blue circles depict the blank fields, and the black sails are the objects of 
the NGC 2000.0 in the region. 



adin, we have provided TESELA with the full capacity and power 
oftheVO. 

In Fig. [7] we show an example of how TESELA visualizes the 
data using Aladin. TESELA sends to Aladin the Tycho-2 stars of 
the region, which are loaded in the first plane. The search area is 
plotted in a second layer with a red circle; another layer depicting 
the blank fields with blue circles is created, and finally, a last layer 
shows the objects of NGC 2000.0 (The Complete New General Cat- 
alogue and Index Catalogue of Nebulae and Star Clusters; [Sinnott| 
1997). Using this Aladin window, users are allowed to load images 



and catalogues both locally and from the VO, having full access to 
the whole universe resident in the VO. This can be very helpful to 
determine the potential influence of relatively bright nebulae and 
extragalactic sources in the sought regions. Note also that, for ob- 
vious reasons, Solar System objects have been not considered, and 
they must be taken into account in order to make use of blank field 
regions close to the Ecliptic at a given date. 

So far, the current version of TESELA allows users to access 
a collection of blank fields obtained from the optical Tycho-2 cata- 
logue. But Tesela has been conceived as a dynamic tool, which will 
be improved in the future with both deeper optical catalogues and 
catalogues in others wavelength ranges. Any future change in the 
tool will be properly documented in its web site. 
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APPENDIX A: HANDLING BOUNDARY TRIANGLES 

Although tessellating small sky subregions is a good approach to 
the problem of deriving blank fields when dealing with very large 
stellar catalogues, the discussion of Section [2.2. 1 1 revealed that the 
Delaunay triangulation of such subregions leads to the situation dis- 
played in Fig. [3J3, where many circumcircles associated to bound- 
ary triangles (shown in magenta colour in that figure) clearly ex- 
pand beyond the limit of the considered subregion. Although in 
this case those circumcircles can be easily discarded, it is possible 
to derive, starting from the coordinates of the corresponding bound- 
ary triangles, new circles devoided of stars that remain within the 
boundary of the sky subregion. This appendix describes in more de- 
tail different approaches that can be employed in this circumstance. 
In what follows we are considering the celestial sphere as the 
unit sphere centered at the origin O. The sky subregion where the 
blank fields will be searched for is the surface of a spherical cap, 
which center is given by the unit vector X S k y , and the angular radius 
of such subregion, as measured from O, is max (see Fig. |Al| ). The 
Delaunay triangulation will be computed with all the stars (down 
to a given magnitude) within that sky subregion. Denoting the unit 
vectors poking to these stars as Si (with i — 1, . . . , A^ sta rs, being 
Altars the total number of stars within the subregion), it is obvious 
that the condition that all these stars must satisfy is 



#max > arccos(K sky -Si), Vz = 1, 



,iVst 



(Al) 



We are also assuming that max < tt/2, i.e., the sky subregion area 
under consideration is smaller than a single hemisphere. 



Discovering problematic boundary circumcircles 

The Delaunay triangulation provides the coordinates of the three 
nodes (stars) of all the triangles resulting from the tessellation pro- 
cedure (see Fig. |Al| ). Let denote the unit vectors pointing from O 
toward the three nodes of each triangle, given in counterclockwise 
order, as Qi , Q 2 , Q3, The unit vector C marking the circumcentre 
of any triangle can be easily computed as (see e.g. Renka 1997 ) 



(Q2-Q1) x (Q3-Q1) 
||(Q 2 -Qi)x(Q 3 -Qi) 



(A2) 



Note that since the three nodes define a plane, the two vectors 
(Q2 — Qi) and (Qs — Q\) are contained in that plane, and their 
cross product provides a new vector perpendicular to the plane. Of 
the two possible directions of that vector, the fact that the nodes are 
given in counterclockwise order guarantees that C points toward 
the spherical triangle defined by Qi, Q 2 , and Q3. In addition, the 
modulus calculated in the denominator of Eq. |A2| assures that C 
is normalised. It is not difficult to show that the circumcentre is 
equidistant from all the three nodes, which allows to compute the 
circumradius of the circumcircle from the inner product of C with 
any of the node vectors 



re — arccos(C • Qi), with i — 1, 2 or 3. 



(A3) 



On the other hand, the angular distance between the circumcenter 
and the center of the sky subregion can be determined as 

dc — arccos(C • K^ y ). 

The condition that must be met for any circumcircle to be inscribed 
within that subregion is then 



If the above condition does not hold, the computed circumcircle ex- 
pands beyond the sky subregion and cannot be used as a valid blank 
field. Under this circumstance other options need to be considered, 
as explained in the next steps. 



Step 0: compute the circle inscribed in the triangle 

No matter how close to the boundary of the sky subregion is located 
any triangle computed during the tessellation, the circle inscribed 
in those triangles are always within the subregion. These incircles 
are centered at a particular point, the incenter, which is equidis- 
tant from the three sides of the triangles and that corresponds to 
the intersection of the three bisectors through the vertices of those 
triangles (see Fig. |A2fr ). The unit vector I pointing from O toward 
the incenter of a given triangle can be obtained as 



1 = 



(JVi >2 ~ AT 3 ,i) x (N 2 ,s - A/i^) 
|(JVi, a -JV 8 ,i)x (N 2 ,s-N 1 , 2 ) 



(A5) 



where AT i?<7 - is the unitary vector perpendicular to the plane defined 
by the origin and the nodes at Qi and Qj, 



N t j = 



Qi x Qi 

I y°£i X C^jf 



(A6) 



Defined in this way, and considering that the three triangle nodes 
are given in counterclockwise order, the three normal vectors A/1,2, 
AT 2 ,3 and AT 3?1 are always pointing toward the interior of the trian- 
gles. It is also easy to show that the inner product I • AT i?:7 - provides 
the radius of the incircle as 



77 = arcsin(i" • Nij), 



(Al) 



whose value is the same for (7, j) = (1, 2), (2, 3) or (3, 1). Note 
the use of the arcsin function instead of the arccos in the previous 
expression; since Nij is perpendicular to the plane defined by the 
origin O and the nodes at Qi and Qj, the angular distance 77 be- 
tween I and that plane must be computed as tt/2 minus the angle 
subtended by I and Nij . 

By construction it is obvious that incircles will always be re- 
gions devoided of stars, since they are inscribed in the considered 
triangles, that in turn are inscribed in the circumcircles, which, by 
the empty circumcircle interior property (see Section |2.1| ) do not 
contain any star. For that reason they constitute a first solution to 
the problem with the boundary triangles which circumcircles ex- 
pand beyond the limit of the sky subregion under tessellation. 



Step 1: enlarge the radius toward the triangle nodes 

Once the incircle of a given boundary triangle has been computed, 
it is very likely that the incenter does not coincide with the circum- 
center. In this case it is possible to enlarge the incircle by increasing 
the radius 77 to the closest triangle node as measured from the in- 
center (see Fig. |A2| 3). To facilitate the explanation of the procedure, 
let rename the unit vectors pointing toward the three nodes of each 
triangle as Qi st , Qind, and Q^d, where the subscript indicates the 
ordering of the nodes as a function of their proximity to the incen- 
ter. 

The enlargement of the angular distance connecting I with 
Qist provides a new radius than can be expressed as 



77, enlarged = arCCOs(J • Qi 



(A8) 



d c +r c < 6m 



(A4) 



However, although the the new enlarged circle will only pass 
through one of the nodes and leave outside the other two nodes, it 
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Figure Al. Graphical illustration of the notation employed in Appendix [A] to explain the initial steps when handling boundary triangles. Diagram (a): The 
whole celestial sphere of unit radius and centered at the origin O is displayed in gray. The sky subregion under tessellation, which boundary is shown in red, 
is centered at the point indicated by the unit vector K s ty, and has a radius 6 max (dashed red line) measured along a great circle on the surface of the celestial 
sphere. Only the stars contained within this region are employed to compute the Delaunay triangulation (shown in cyan). Diagram (b): A sample boundary 
triangle formed by the stars located at Qi, Q2 and Q3, given in counterclockwise order, is represented in green. Its corresponding circumcircle, centered at 
C and passing through the three node stars, is displayed in magenta. The radius of this circumcircle re is plotted with a dashed magenta line. The angular 
distance dc between C and -K~ s ky is shown with a dashed blue line. In this example the circumcircle expands beyond the limit of the sky subregion, i.e. the 
condition given in Eq.|A4|does not hold. 



is possible that other neighbouring stars can enter into its interior. 
For that reason, instead of Eq.|A8|it is necessary to use 



ri, enlarged = minimum [arccos(J ■£»)], Vz = 1, 



,iVst, 



(A9) 



which already includes the relation given in Eq. |A8[ since the three 
nodes are part ot the set of stars. 

Before accepting r/, e niarged as a valid radius for a blank field, 
it is required to check that with this new radius the resulting circle 
(centered at the triangle incenter) does not exceed the boundary of 
the sky subregion. This is performed in a similar way as previously 
employed in step 0. First the angular distance between the incenter 
and the center of the sky subregion is determined as 

di — arccos(i" • K^ y ). 

The condition that must be met by the enlarged circle is then 



*>I \ ^ I, enlarged S: 



(A10) 



If the above condition is not satisfied, a revised version of the radius 
can be computed as 



f I , enlarged, revised 



d T . 



(All) 



Note that, so far, the valid blank field is still centered around 
the triangle incenter. However, by removing this constraint, it is 
possible to obtain an even larger valid blank field by allowing its 
center to move from the incenter toward the directions of the other 
two nodes following the great circles connecting these points. 

Next we describe the procedure of shifting the centre of the 
blank field along the line connecting the incenter with the node 
corresponding to Q2nd- The resulting new center V (see Fig. |A2| :) 
will be a point equidistant from Q\ st and Q2nd- The situation with 



the third node will be identical interchanging Q2nd by Q3rd in the 
following description (see Fig. |A2fl ). 

The first step consists in computing the midpoint in the tri- 
angle side connecting Q\ st and Q2nd- The unitary vector pointing 
toward this point is given by 



M 



Qist + Qi 

IQist + Qi 



(A12) 



The sought new center V can be obtained as the intersection of the 
great circle defined by two vectors M and C with the great circle 
defined by the two vectors I and Q2nd- This intersection can be 
easily computed as 



V ■ 



(CxM)x(Ix Qi 
\{C x M)x (Jx Q 2 



(A13) 



where the ± sign indicates two initial solutions from which the 
closest vector to i" must be selected. The angular distance from the 
new center V to either Q\ st and Q2nd is now 



r v 



arccos(V • Q no de), with node=lst or 2nd. 



(A14) 



Next, it is necessary to check that no new stars (different from the 
ones at the triangle nodes) have entered within the new circle. This 
can be computed as 

?V,revised = maximum {?v, minimum [arccos(V • Si)]} , (A15) 

where i = 1, . . . ,7V stars . 

Note that it must also be checked that with this new radius 
prevised the circle centered at V is still within the sky subregion. 
In a similar same as we proceeded in Step 2, the distance from V 
to Xsky is determined as 
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Figure A2. Different strategies described in Appendix|A|to obtain a valid blank field in the neighborhood of a boundry triangle. These diagrams represent an 
enlargement of the example represented in Fig. |Al| where, for clarity, some graphical elements represented in that figure have been removed. 

— Diagram (a): Computation of the incircle, the circle inscribed in the spherical triangle. It is centered at J and is displayed in blue. The distance from i" to 
each of the sides of the boundary triangle is the same, r j . 

— Diagram (b): Enlargement of the initial incircle until reaching the closest star in the boundary triangle. Note that the three stars at the vertices of the 
spherical triangle have been renamed as Q\ s t, Qind and Q?> x & to indicate their order of proximity to the incenter. Thus the resulting new radius, 77 , enlarged* is 
the angular distance between the incircle J and Qi s t- 

— Diagram (c): Displacement of the center of the blank field along the great circle connecting the incenter i" with the second nearest triangle node Q2nd- The 
new center V is equidistant to the nodes Qi st and Q2nd- 

— Diagram (d): Displacement of the center of the blank field along the great circle connecting the incenter I with the furthest triangle node Q3 r d- The new 
center V is equidistant to the nodes Qi st and Q3 T d- 

— Diagram (e): Displacement of the center of the blank field along the great circle connecting the incenter i" and the circumcenter C. The distance rw from 
the new center W to the nearest triangle star is the same than the distance to the boundary of the sky subregion, indicated by F. 

— Diagram (f): Random search of a valid blank field around the spherical triangle defined by Qi, Q2 and Q3. The black dots are random points in the 
intersection of the circumcircle and the sky subregion. Each of these points is checked as the center of a potential blank field. The blue circle, with center T 
and radius rj>, indicates the largest of such circles that satisfy the conditions for a valid blank field. 



dv — arccos(V • i*T S k y ), (A 16) 

and then the condition that must be satified can be written as 



dy + ?V,revised < #m 



(A17) 



If the previous condition does not hold, the radius of the blank field 
can be redefined as 



^V,revised2 — "max &V • 



(A18) 



The resulting blank field centered at V will be a better (larger) 
blank field as far as revised (or 7v jrev ised2) is larger than 77, enlarged- 



Step 2: shift the center between the incenter and the 
circumcenter 

Independently of the success of the previous step, another possi- 
bility worth being explored is shifting the center of the blank field 
along the great circle connecting the incenter with the circumcen- 
ter, with the constraint that the new circle does not cross the sky 
subregion boundary (see Fig.|A2fr). 



For this computation the bisection method can easily be em- 
ployed. Since the solution for the new center will be a point in the 
great circle connecting i" with C, it is possible to use two auxiliary 
vectors L\ and L2 which are initialized as Li —I and L2 = C. 
Then the process enters into an iterative procedure in which the 
following steps are performed: 
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(i) The midpoint on the great circle connecting L\ and L 2 is 
determined as 



-t'mid — 






(ii) The minimum angular distance from the previous point to 
any of the triangle nodes is evaluated using 

g?l-q = minimum [arccos(L m id • Q*)], with i — 1, 2 or 3. 

(iii) The angular distance from Lmid to the border of the sky 
subregion is also determined as 



^L-Sky 



arccosfXn 



•K t 



sky ) 



(iv) If the two previous angular distances are equal within a 
given tolerance (for our purposes |c2l-q — c^-sky! < 0.01 arcsec 
is enough), the iterative procedure halts and the final result is de- 
fined slsW — Lrmd. Otherwise one of the auxiliary vectors must be 
redefined according to 



> 0, then L 2 = L m [ d 
< 0, then Li = L mid 



if (cL-q - d L -sk y ) is • 

and the process is iterated by repeating steps (i)-(iv). 



At the end of the iterative process a new center W has 
been computed (see Fig. |A2fr ), which is equidistant to the near- 
est triangle node and to the boundary of the sky subregion, i.e., 

rw — i-Sky = <^L-Q- 

Finally it must be checked that no new stars have entered into 
the new circle. Similarly to what was done in previous steps, the 
new radius can be refined by using 

rw,revised = maximum{rv^ , minimum [arccos(W • Si)]} (A 19) 

where i = 1, . . . , N^s- 

It is important to highlight that new radius r^revised will not 
necessarily be larger than the radii previously derived in Steps 2 
and 3 (in fact, the example illustrating this appendix shows that the 



blank field in Fig. |A2 
the blank field in Fig. 



I centered at V with radius rv is larger than 
A2fe centered at W with radius rw)- 



Step 3: random search in the intersection between the 
circumcircle and the sky subregion 

A final, but also effective, brute force method consists in using a 
random search of a valid blank field center within the intersec- 
tion region between the circumcircle and the sky subregion (see 
Fig. |A2f ). For this purpose, random points are generated on that 
region, imposing that any small area has to contain, on average, the 
same number of points (see e.g. Weisstein 201 1 ). For each random 
point T the valid blank field radius is evaluated as 



r T = minimum [arccos(T • Si)], \/i = 1, 



,N st 



(A20) 



It is also necessary to check that the resulting blank field remains 
inscribed in the sky subregion. For that purpose the distance from 
the center to the border of the subregion is determined using 



dr — arccos(T • X" S ky), 

and the corresponding condition to be verified is 

d T +TT < ma x. 



(A21) 



(A22) 
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Number of random points in step 3 

Figure A3. Box-and- whisker plot showing the percentage of success of the 
different steps described in Appendix [A] The bottom and top of each box 
show the 25th and 75th percentiles, whereas the thick horizontal line marks 
the median. The whiskers (vertical dashed lines) show 1.5 times the in- 
terquartile range. Since each colour represents a fixed number of random 
points employed in the search for valid a blank field region in step 3, the 
sum of the percentages in the three steps for a fixed colour is 100. 



Obviously, the solution obtained with this method must be 
compared with the solutions derived in the previous steps in order 
to decide which one provides the largest blank field in the neigh- 
borhood of the considered boundary triangle. 



Successfulness of the different steps 

The three steps just described produce different possible valid blank 
fields associated to each boundary triangle. The adopted valid blank 
field in each case will be the largest among those alternatives. 

It is interesting to investigate the successfullness of the three 
steps providing the largest blank field. As a first-order estimation of 
these numbers, we have measured the percentage of success of each 
step when tessellating random stellar fields. Since the success ratio 
for step 3 depends on the number of random points drawn within 
the intersection region between the circumcircle and the sky sub- 
region, we have computed those percentages using different num- 
bers of random points. The results are displayed in Fig. |A3| Each 
colour in this figure indicates a fixed number of random points em- 
ployed in step 3. When no random points are used (represented in 
red colour), the median percentage of success of steps 1 and 2 are 
23% and 77%, respectively. Not surprisingly these values decrease 
as the number of random points in step 3 increases, reaching a sta- 
ble situation as far as the number of random points exceeds a few 
hundred. In particular, for 10 4 random points the median percent- 
ages of success (represented in magenta colour) are 10%, 24% and 
66% for steps 1, 2 and 3, respectively. 

The previous results have shown that although step 3 will be, 
in general, the resposible for providing the valid blank field, the 
contribution of the other two steps is not negligible. 



If this is not the case, the value of the radius can be refined using 



^T,refined = ' 



dr. 



(A23) 



